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Metadynamics is an established sampling method aimed at reconstructing the free-energy surface 
relative to a set of appropriately chosen collective variables. In standard metadynamics the free- 
energy surface is filled by the addition of Gaussian potentials of pre-assigned and typically diagonal 
covariance. Asymptotically the free-energy surface is proportional to the bias deposited. Here 
we consider the possibility of using Gaussians whose variance is adjusted on the fly to the local 
properties of the free-energy surface. We suggest two different prescriptions: one is based on the 
local diffusivity and the other on the local geometrical properties. We further examine the problem 
of extracting the free-energy surface when using adaptive Gaussians. We show that the standard 
relation between the bias and the free energy does not hold. In the limit of narrow Gaussians 
an explicit correction can be evaluated. In the general case we propose to use instead a relation 
between bias and free energy borrowed from umbrella sampling. This relation holds for all kinds of 
incrementally deposited bias. We illustrate on the case of alanine dipeptide the advantage of using 
adaptive Gaussians in conjunction with the new free-energy estimator both in terms of accuracy 
and speed of convergence. 



I. INTRODUCTION 

The problem of sampling complex energy surfaces, 
characterized by metastable states separated by large en- 
ergy barriers, has recently received considerable atten- 
tion. A list that is by no means exhaustive of the pos- 
sible remedies suggested includes transition path sam- 
pling PP, umbrella sampling [3], local elevation [3], 
Wang-Landau [5] , adaptive biasing force [5] , metadynam- 
ics [7HTU] and self healing umbrella sampling [TT] . 

We focus here on metadynamics (MetaD) [5j. In 
this approach one starts by identifying a set of appropri- 
ately chosen collective variables s which are function of 
the microscopic variables q. In order to accelerate sam- 
pling, a bias potential is dynamically added during the 
simulation. The bias has the effect of helping the system 
to overcome large free-energy barriers so as to accelerate 
sampling. Asymptotically, the negative of the bias pro- 
vides an estimate of the free energy F(s) associated with 
the collective variable s. The estimate has been demon- 
strated to be free of systematic errors, if the CVs are 
properly chosen [S]. 

Metadynamics has been successfully used in several 
different contexts, for a recent review see for instance 
Ref. [ID]. In most cases the bias is constructed by peri- 
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odically adding a repulsive Gaussian potential [fj which 
is a function of the s. This implies defining the height and 
width of the Gaussians. The first is related to the energy 
deposition rate while the latter is taken small enough to 
resolve the free-energy surface features. Prescriptions on 
how to choose these parameters have been given and the 
dependence of the statistical error on their choice has 
been discussed [8] . Still one is not guaranteed a priori of 
making the best choice. 

More recently, a new flavor of metadynamics has been 
introduced which goes under the name of well-tempered 
metadynamics (WTMetaD) .[12] In WTMetaD the speed 
at which the bias is added decreases during the simula- 
tion. WTMetaD maintains the property that the asymp- 
totic bias is related to F(s) by a simple relation but, 
at variance with standard metadynamics, the final free- 
energy estimate converges to a definite limit. 

A useful property of WTMetaD is that the dependence 
of the final result on the speed with which the bias grows 
is smaller than in standard metadynamics |12j . reducing 
the impact of an improper choice of this parameter. Here 
we want to move further towards making the method 
even more efficient and robust with respect to parame- 
ters choice. A natural step in this direction is to add the 
possibility of adapting the Gaussian width to the local 
free energy so as to speedup sampling. For instance a 
free-energy surface may present minima with rather dif- 
ferent curvature and a Gaussian width optimal for one 
might be not appropriate for another (see Supplemen- 
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tary Information, Section [Dj) . This problem is amplified 
in the case of CVs that are are highly non-linear functions 
of the atomic coordinates. 

In the past this need was recognized. Therefore ad 
hoc solutions for speeding up free-energy surface explo- 
ration were developed [13] and some form of metady- 
namics that uses adaptive Gaussians was implemented 
in publicly available codes [T3]. Very recently Tribello et 
al [15] also used adaptive Gaussians in a related but dis- 
tinct context. However, the effect of these choices on the 
free-energy reconstruction has never been systematically 
investigated. 

Another common assumption in multidimensional 
MetaD simulations is to adopt Gaussian functions whose 
axes are aligned to the chosen CVs, i.e. with a diagonal 
covariance matrix. This is clearly a simplification which 
reduces the number of input parameters, but an opti- 
mal choice might require Gaussian axes which are not 
aligned to the CVs. This was pointed out in Ref. [T3] 
and successively has been used for crystal-structure pre- 
diction [23 [T7] . 

In this paper, we discuss the possibility of performing 
metadynamics simulations using multivariate Gaussian 
potentials with a full covariance matrix computed on the 
fly. We first introduce two possible different schemes for 
choosing the shape and width of the Gaussian poten- 
tials. One is dynamical and based on the mean square 
displacement of the CVs in a predetermined time inter- 
val. The other is geometrical and based on the local 
mean square displacement of the CVs due to a hypo- 
thetical change in the microscopic coordinates. Subse- 
quently we observe that the adoption of either of these 
two position-dependent covariances requires modifying 
the free-energy estimator in order to reconstruct accu- 
rately the free-energy landscape. Furthermore, we show 
that this new scheme applies also when the added Gaus- 
sian potentials or some other chosen form of potentials 
are coarser than the underlying free energy. The algo- 
rithms and their efficiency are numerically tested in the 
standard case of alanine dipeptide in vacuum. Details 
of the simulations and error calculations arc discussed in 
the Appendix. 



II. METHODS 

A. Metadynamics 

In MetaD, a time-dependent bias potential V(s, t) is 
used to discourage the system from visiting already ex- 
plored regions in the CVs space [7J [T5] . At the beginning 
of the simulation the bias is everywhere equal to zero. 
Subsequently it is evolved according to the following ex- 
pression: 

where g(s, s(t)) is a short-ranged kernel function. 



Here AT is an energy which can be used to tune the 
region of free energy explored, and u> is the initial filling 
rate. It has been shown in Ref. [T^] that for large times 
the bias is related to the free energy by 

F(s) = - hm *££L\y(s,t)-C(t)] . (2) 

where C(t) does not depend on s. We notice here that 
even if C{t) eventually diverges, its impact on the result is 
completely irrelevant as it disappears in any free-energy 
difference calculation. The limiting behavior of C(t) is 
discussed in Section 

These two expressions covers both the case of 
WTMetaD[12] (where AT is finite and is an additional 
tuning parameter) and non-WTMetaD[7J, which is sim- 
ply recovered in the limit AT — > oo and corresponds to 
the case where the filling rate is kept constant during the 
simulation. The properties of Equation [TJ were exten- 
sively discussed in Ref. [15] . Equation [2] is of practical 
use since it offers an estimate of the free energy up to 
an arbitrary constant if one assumes its validity also at 
finite time 

FM = - AT ^[V(s,t)-C{t)] . (3) 

In most practical applications, the kernel g(s, s(t)) of 
Equation [TJ has been replaced by a smooth function of the 
CVs, typically a product of one dimensional Gaussians. 
We choose here a generic multivariate Gaussian function: 

GX P [~\ X> - M*)K 2 [ Sj - sj(t)]j . (4) 

Here, cry determines the shape of the Gaussian poten- 
tial. In conventional MetaD simulations, o~ij have almost 
always been assumed diagonal with only few exceptions 
[TBI fTBTtlT] . We stress once again that in all these cases 
MetaD was not employed to estimate the underlying free 
energy, but only to speedup the exploration of configu- 
ration space and no attention has been paid to the effect 
of these protocols on the free-energy reconstruction. 

Passing from a diagonal <xy to a non-diagonal one in- 
creases the number of parameters that have to be chosen 
when setting up the simulation, so that it would be con- 
venient to give a prescription that simplifies this choice. 
In standard MetaD the diagonal elements of cr.y are ob- 
tained by computing the CV standard deviation during 
a short preliminary simulation. A simple extension is to 
consider also the correlation between different variables 

o-?. oc (AsiAsj) . (5) 

This choice is completely invariant with respect to an 
arbitrary linear transformation in the CV space and re- 
duces to a diagonal try by a suitable transformation [16j . 
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In addition it allows easy mixing of CVs of different units 
and nature. However an optimal (Tij over the entire CVs 
domain requires that is made position dependent. Thus 
we shall still use Equation [5] to estimate the Gaussian 
covariance but, in the same spirit of Ref. [13] , we shall 
give a time dependent estimation of (As^Asj) so as to 
reflect the local properties of the free energy. We finally 
note that the use of Equation [5] still preserves the prop- 
erty of invariance relative to linear combination of the s 
variables. However, the covariance matrix cannot be re- 
duced to a diagonal form everywhere by means of a single 
linear transformation of the CVs since the matrices ay at 
different times or positions cannot be expected in general 
to commute. 



B. Dynamically-adapted Gaussian 

To define a time-dependent adaptive covariance at time 
t we compute the average value of the CVs and all the el- 
ements of the covariance matrix from the last part of the 
trajectory. The center of the Gaussian is placed at the 
computed average value. To select the segment of tra- 
jectory over which we perform the average we found con- 
venient to introduce an exponential weighting function 
with characteristic decay time td, such that the Gaus- 
sians' centers Sj(t) and their covariances at time t are 
given by: 



8i(t) 



1 



T D Jo 



dt' Si (t') e -(*-*'V^ 



(6) 



and 



dt%{t')-s i (t')][s j (t')-s j (t')} e -(*-'')/r D . 

(7) 

With this choice the Gaussian location and covariance 
change smoothly and can be very easily evaluated. To 
this effect we take the time derivative of Equation [6] and 
Equation [7] 



i(t) = 



i(t)-3j(t) 

td 



(8) 



Within this scheme td determines the time window 
which is used to estimate the CV fluctuations and thus 
to choose the Gaussian width. It is instructive to see 
how the latter depends on the dynamical properties of 
the system at least in the simplified case of Langevin dy- 
namics (see example in the Supplementary Information, 
Section [D]). In this case, two regimes can be identified, 
for short and long value of the simulation time. At the 
beginning of the simulation the dynamics is still stuck in 
the metastable minima. If td is larger than the typical 
autocorrelation time of the CV the Gaussian shape will 
be equal to the shape of the corresponding free-energy 
minimum, thus providing an optimal filling. At the end 
of the simulation, when the barriers have been smoothed 
out by the adaptive bias, the dynamics is close to a free 
diffusion, with a diffusion coefficient which is possibly po- 
sition dependent. Simple dimensional considerations can 
be used to show that the computed afj matrix becomes 
proportional to the position-dependent diffusion tensor 
Dij(s(t)), which has a relevant role in describing many 
important phenomena. [18j 



C. Geometry-adapted Gaussians 

We shall now discuss an alternative protocol for choos- 
ing the a matrices. Our starting point will be again 
Equation [5] For small displacements of the microscopic 
variables q the associated change in each CV can be lin- 
early approximated as 



dsj 
dq a 



Aq a 



(10) 



Then in this approximation it can be easily seen that 
the CVs covariance is linearly related to that of the 
atomic displacements (Aq a Aqp). If we assume these to 
be Gaussian distributed with standard deviation oq we 
have (Ag Q Ag^) = & a po-Q leading to 



(As,- As,- 



x - dsi dsj 
„ dq a dq a ' 



(11) 



We use this expression to define the shape of the Gaussian 
covariance via the Gram matrix: 



/(*) = 



»(*) -«<(*)] 



%(*)]- 



td 



(9) 



and then consider Sj(t) and erf At) as additional variables 
to be evolved together with the system dynamics. Inte- 
grating Equation [8] and Equation[9]with the initial condi- 
tions Sj(0) = Si(0) and cr^(0) = the values of Si(t) and 
afj(t) as defined by Equation |6| and Equation [7] are recov- 
ered. We name this scheme aynamically- adapted (DA). 
Using partial time averages to determine u 2 {i) is some- 
what natural and it has the practical benefit that only 
one parameter td determines the whole covariance ma- 
trix. 



M) = °~G 



2 v OSi dsj 



dq a dq a 



(12) 



With this choice, cr 2 (q) depends explicitly on the micro- 
scopic variables q. We shall refer to this case as geometry- 
adapted (GA). As in the DA case only one parameter suf- 
fices to determine the whole covariance matrix. We note 
that in MetaD the Gaussian spread plays a role similar 
to that of the histogram bin size for other methods and 
determines which configurations are considered as equiv- 
alent. Thus here we are assuming as equivalent all the 
microscopic configurations whose a root square distance 
of the atoms involved in the CVs is within o~q making 



4 



possible for the choice of uq to be guided by physical con- 
siderations. For instance, if one were to choose a different 
set of CVs still the choice of sigma would be determined 
by the typical atomic displacements. 

Similarly to the DA scheme, the GA one is invariant 
with respect to linear transformations of the CVs. 



III. FREE-ENERGY ESTIMATION 

One of the most interesting features of metadynamics 
is the link it establishes between the free-energy surface 
and the bias both in its classical and well-tempered ver- 
sions. Having changed the protocol of Gaussians depo- 
sition it is crucial to establish such a link in the case of 
Gaussians of variable covariance. 



A. Small-width limit 

We consider first the case in which the Gaussians have 
a size smaller than the relevant features in the free-energy 
surface, as normally done in MetaD. This was the as- 
sumption of Ref. [12] where one went as far as to assim- 
ilate the Gaussians to 8 functions. 

Within our schemes this limit can be achieved by 
choosing tjj or <jq small enough. In this regime, for the 
DA scheme, one can also neglect the difference between 
the actual position s and the center of the deposited 
Gaussian s. Therefore the only adjustment necessary is 
to take into account the fact that the change in covariance 
induces a change in Gaussians' volume. Thus in Equa- 
tion [I] we replace the kernel g(s,s(t)) with a Dirac delta 
with the same normalization, \J (27r) d det aS(s — s(t)), 
where a is either a function of the trajectory in CVs' 
space (DA scheme) or of the microscopic coordinates at 
a specific time (GA scheme) and d is the number of CVs: 

V(s, t) = we -V{ S (t),i)/AT^ 2?r ^ det aS ( s _ (13) 

As in Ref. [T2], we notice that for large times 
the probability distribution becomes P(s,i)ds oc 
exp ( — - F ( s ) H ^ / ( s ' t ) j anc | one has 



This means that the bias is no longer proportional to 
F(s) as in Equation [2] but rather to 



V(s, t) = we- y(s,t)/AT ^/(27r) d (det a) s P{s, t) 



-V(s,()/AT 



(det a) 



-[F(s) + V(s,t)]/T 



(14) 



where the average (dettr) s is taken in the canonical en- 
semble at a fixed value of the Gaussian center. Now, 
setting V — C(t) in the last equation, where C{t) is con- 
stant with respect to the CVs, we find that the asymp- 
totic solution for V(s,t) is 

Km V(s, t) = - [F(e)-T ln(det a) s ] + C(t). 

(15) 



G(s) = F(s) -Tln(deto-) 



(16) 



In the case of GA scheme of Section II C G(s) turns 
out to be the gauge invariant free-energy discussed in 
Refs. [19-21 which is not changed by an arbitrary non- 
linear monotonic transformation of the CVs. Instead in 
the DA case of Section |IIB| when td is chosen so that a 
is proportional to the diffusion matrix, our bias poten- 
tial is closely connected with that used in flux-tempered 
metadynamics [22 which is also gauge invariant and has 
been designed to optimize the round trip time. 

While gauge invariance is esthetically pleasing it does 
not bring particular advantages since for physical appli- 
cations it is F(s) that is needed. This can be simply 
obtained by rewriting Equation |15| as 



lim A ^ T [V(s, t) - C{t)] + T ln(dct a) , 



F(s) 

(17) 

which is the appropriate generalization of Equation [2] to 
the case of small adaptive Gaussians. 



B. Free energies from reweighting 

As discussed in the previous sections, for non-standard 
biasing protocols the relationship between the asymp- 
totic bias and the underlying free-energy landscape is 
not known a priori. Only in the small width limit it 
is possible to estimate explicitly the correction of Equa- 
tion 1171 However it would be nice to have an estimator 
that is valid for large Gaussians. This is provided by the 
relation: 

F N (s, t) = -T In N(s,t) -V(s,t) + T log J ds'N(s',t) 

(18) 

where N(s,t) is the accumulated histogram of the vari- 
able s up to time t. If the bias is time-independent 
this relation is strictly true and is normally used in um- 
brella sampling. In this present context it is valid only 
for large times when V(s,t) has converged. More pre- 
cisely it is only necessary that the rate of bias deposition 
goes to zero sufficiently rapidly. Therefore it will be valid 
whether the Gaussian width is small or large and whether 
the width is determined by diffusion or geometry. 

This relation was also used in deriving a free-energy 
estimator for WTMetaD [T^] where this expression was 
manipulated in a manner similar to Section [ill A| leading 
to the estimator in Equation^ However, F/v(s, t) is more 
generally valid and gives a correct estimate of free energy 
even when adaptive Gaussians are used. In the following 
we will show how this new estimator can be of use in a 
practical case. 

The case of non-WT metadynamics (i.e. for AT — > 
oo ) where the bias oscillates in the long time limit needs 
separate consideration and will be discussed elsewhere. 




FIG. 1: Molecular sketch of alanine dipeptide with the Ra- 
machandran dihedral angles $ and This graphics was pro- 
duced with VMDG3. 



IV. EXAMPLES 



FIG. 2: The estimate of the average error as a function of 
simulation time for alanine dipeptide in vacuum. The average 
error was calculated through a set of 100 runs and av erag ed 

for 
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every 120 ps. The effect of adopting Fn of Equation 
two different choices of a is compared to the choice of F of 
Equation [3] It is possible to appreciate that Equation [18] al- 
ways delivers results that are comparable to the Equation [3] 
or better, when large Gaussians are adopted. In the inset 
it is shown the time dependence of the error that decreases 
linearly with 1 / y/t predicting an error that goes to zero for in- 
finite time. On the contrary, for large Gaussians, the standard 
estimate is not able to resolve the features of the free energy 
resulting in a significant residual error even for infinite time. 



Narrow and wide Gaussians 



Before discussing the performance of the adaptive 
Gaussians we shall first validate the estimator of Equa- 



tion 18 when using standard WTMetaD and compare 
two very different choices of Gaussian width. In this way 
we show that it is advantageous to use Equation 1 1 8| even 
when using Gaussians with fixed covariance. 

To this effect we performed a number of WTMetaD 
simulations of alanine dipeptide in vacuum (see Figure fl]) 
using as collective variables the two Ramachandran |z4j 
dihedral angles $ and ^ . The initial Gaussian height 
was chosen to be 0.287 kcal mol -1 and Gaussian poten- 
tials were deposited every 120 fs, corresponding to an 
energy deposition rate of 2.39 10 -3 kcal mol -1 fs -1 . The 
AT parameter for WTMetaD was set to 1200 K. Other 
technical details can be found in Section FBI 

We first performed a long calculation to obtain a refer- 
ence free-energy landscape. This was done in two steps: 
we performed a 5 ns long WTMetaD run with a Gaus- 
sian width cr= 0.35 rad for both $ and ^ and the bias 
thus accumulated was then kept constant in a long (1/is) 
biased simulation. Then, by using Equation |18| which is 
exact for a static bias, we obtained the reference free en- 
ergy. The error of the reference landscape was estimated 
by comparing the free energy derived from the histogram 
from the first and the second half of the simulation and 
was obtained using the procedure reported in Section |C| 
The error is approximately 0.01 kcal mol -1 , which is neg- 
ligible as compared to the error of more than 0.05 kcal 
mol -1 made in the test runs described below. We did 
not bin the s to evaluate the histogram iV(s,t) in Equa- 
tion [18] but instead we used Gaussian functions of width 



0.01 rad. This led to a smooth N(s, t) and a smooth free 
energy as shown in Figure [4]^. 

We then performed two different sets of simulations 
using the same a for both dihedral angles. In one case a 
was smaller than the free-energy features (a= 0.35 rad), 
in the other larger (cr= 0.7 rad). Each set of simulations 
consisted of 100 WTMetaD runs of 20 ns each. The error 
relative to the reference free energy committed using the 
standard estimator of Equation [3] and the new one in 
Equation 1 18| was then compared. 

From Figure [2] it can be seen that in all cases the new 
estimator produces consistent results and the expected 
asymptotic 1/yt behavior. The old one works well for 
small Gaussians but in the case of the large ones is af- 
fected by a systematic error, due to the fact that it is 
not able to resolve the smaller features of the free energy 
landscape. 



B. Dynamically-adapted Gaussians 

Having assessed the usefulness of Equation [18| we turn 
to the evaluation of the performance of Gaussians of vari- 
able covariance using the DA scheme. This requires defin- 
ing the parameter T£> that is a measure of the time re- 
quired to sample two bins that we consider as different. 
In order to choose tjj in a manner which allows a fair 
comparison of constant and variable Gaussians covari- 
ance we performed a preliminary 1 ns of standard WT- 
MetaD with Gaussians of a constant er=0.35 rad width. 
We then tested different td values and we choose the 
one that was able to fill on average the same volume 
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change in volume of the Gaussians, which measures the 
error made by using the old estimator. It can be seen 
that in this case the maximum error is as small as 0.6 
kcal mol -1 . However this is by no means generally the 
case as we shall sec in the next section. 



D. Geometry-adapted Gaussians: double MSD 



FIG. 3: Calculation of the average error for both GA (de- 
noted with o q ) and DA (denoted with td) schemes for ala- 
nine dipeptide as a function of time. The average error was 
calculated through a set of 100 runs and averaged every 120 
ps. Two free-energy estimates are employed: the one of Equa- 
tion 
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denoted with Fn and the one of Equation [3] denoted 
with F. The performance of Equation 18 is always superior to 
that of Equation [3] which is dramatically incorrect when DA 
or GA are employed. In the inset it can be appreciated that 
the error scales linearly with as expected in WTMetaD. 

The use of Equation [18] leads to the correct asymptotic be- 
havior in contrast to Equation [3] 



deposited in the standard run. In this way we ensured 
that the system was subject to a comparable filling rate. 
This gave the value td= 300 fs. In analogy to the pre- 
vious protocol, we performed 100 runs of 20 ns each and 
calculated the average error made in the region of free- 
energy surface which is within lowest 5 kcal mol -1 from 
the minimum of the reference surface. 

The results are displayed in Figure [3j where the failure 
of the old estimator of Equation [3] is evident. Instead 
the new estimator gives correct results even when the 
covariance of the Gaussians is let to vary. 



C. Geometry-adapted Gaussians: Ramachandran 

plot 

We now turn to a test of the GA scheme for the 
same system using an identical protocol. The choice 
of a comparable ac was performed as before imposing 
the condition of the same volume filling rate with re- 
spect to a preliminary standard WTmetaD. This led to 
a value oq =0.2 A. The results are again displayed in 
Figure [i) The free energy calculated through F of Equa- 
tion [3] is expected to converge to the gauge invariant 
free-energy landscape, which is different from the stan- 
dard one. Thus, even for long times, there is a system- 
atic residual error. On the other hand, the new estima- 
tor does give the correct result, with a convergence rate 
which is very similar to that obtained with the standard 
calculation reported in Figure [2] 

In Figure [3J3, we give a pictorial view of the shape of 
the Gaussian potentials obtained by using Equation |12| 
We also report in Figure |4p the correction due to the 



We examine now the role that a different choice of col- 
lective variables might have on the error made by using 
the old estimator. While for the angle $ and 'J the error 
was small we expect it to be larger when dealing with 
CVs that are strongly nonlinear functions of the atomic 
coordinates. One such example are for instance CVs 
that measure the mean square deviation (MSD) from a 
reference structure or coordination functions that have 
sigmoidal dependence on the distance of two atoms or 
groups of atoms. 

Thus we studied once again alanine dipeptide in vac- 
uum and used as CVs the MSD from the two conformers 
C7e q and C7 ax rather than the two torsional angles. The 
MSDs were calculated through optimal alignment by us- 
ing Kearsley's algorithm and only the heavy atoms 
were considered in the metrics. In order to obtain refer- 
ence values as accurate as possible rather than reweight- 
ing the configurations previously generated, we repeated 
the simulation with the same protocol as before. The 
result of this run which also lasted Ifis is shown in Fig- 
ure [5K where it can be seen that in these new variables 
the free energy surface has two minima one very narrow 
and the other wider. 

By comparing Figure [5)3 and Figure [4^3 it is evident 
that the double MSD space induces a much larger change 
in shapes and volumes of the deposited Gaussian poten- 
tials. Correspondingly the error made using the old esti- 
mator of Equation [3] becomes larger (see Figure[5]C). The 
benefit of using variable Gaussians becomes very appar- 
ent since we do not have to choose small Gaussians to 
resolve the narrow minima paying the price of a slow 
convergence, nor do we have to use larger Gaussians and 
sacrifice accuracy. This is exemplified in Figure [6] where 
we show the typical convergence behavior obtained by 



using the GA scheme of Equation 12 and compare it to 
fixed-cr runs, where a was chosen to be appropriate either 
to the narrow minimum (o=0.01 A 2 ) or to the larger one 
(cr=0.3 A 2 ). Using the new estimator all three calculation 
appear to converge to the same limit, but with different 
rates. For narrow Gaussians (ct=0.01 A 2 ) the calculation 
converges to the right limit albeit rather slowly. Similarly 
for very large Gaussians (o=0.3 A 2 ) the convergence was 
very slow for reasons similar to those discussed earlier 
(see Section IV A | . When the GA scheme of Section II C 
is chosen the asymptotic limit is reached faster and with 
a much smaller error. 
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FIG. 4: Panel A: free-energy landscape on the two Ramachandran dihedral angles <E> and ^. The dotted region represents the 
portion used for error calculations which lies within 5 kcal mol -1 from the lowest free-energy point. Panel B: representative 
sketch of the shape of the Gaussians produced on the Ramachandran plot for alanine dipeptide obtained by using the GA 
scheme of Section |II C| The size of each ellipse is scaled so to reflect the average size of the Gaussian potentials placed in each 
specific point. Panel C: free-energy contribution coming from the change in volume of the Gaussians induced by the adoption 
of the GA scheme. 
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FIG. 5: Panel A: free-energy landscape on the MSD from the two metastable basins C7 oq and C7 ax . The dotted region 
represents the portion used for error calculations which lies within 5 kcal moP 1 from the lowest free energy point. Panel B: 
representative sketch of the shape of the Gaussians produced by using two MSD as CVs. The Gaussian widths were obtained 
by the GA scheme of Section |II C| The size of each ellipse is scaled so to reflect the average size of the Gaussian potentials 
placed in that specific point. Panel C: free-energy contribution from the change in volume of the Gaussians induced by the 
adoption of the GA scheme. The blank region in the center is due to the lack of sampling for this simulation time. 



V. DISCUSSION AND CONCLUSIONS 

In summary, we have investigated in detail the pos- 
sibility of performing metadynamics simulations where 
the repulsive Gaussian potentials have a width chosen 
on the fly and are not aligned with respect to the CVs. 
We have shown that using adaptive widths can lead 
to artifacts in the estimation of the free-energy land- 



scape, which can be recovered by using a suitable es- 
timator. Moreover, we have discussed two indepen- 
dent recipes to adapt the Gaussian shape, one based 
on the time evolution of the collective variables and an- 
other based on the intrinsic metrics of the microscopic 
coordinates. Both methods can be implemented using 
the same ingredients as standard metadynamics calcula- 
tions, namely collective-variable values and first deriva- 




FIG. 6: WTMetaD for the double MSD case for alanine dipep- 
tide in vacuum. For each choice of a, 100 runs of 4 ns each 
were performed and the average error respect to the refer- 
ence free energy was calculated as function of time. Two 
different choices of a are made and compared with the re- 
sult obtained using the GA scheme with Gaussian width oq. 
The GA scheme is always superior and converges faster to the 
correct value. 



tives with respect to atomic positions. They are effective 
and have been shown to produce unbiased results on a 
standard benchmark. Additionally they can remarkably 
improve the filling speed especially when collective vari- 
ables which are highly non-linear functions for the atomic 
positions are used. This is particularly helpful for MSD 
and MSD-based variables or contact map-based vari- 
ables which are widely employed for systems of biological 
and condensed matter interest. For such variables, the 
usual choice of the Gaussian width fixed once and for 
all to the value of the CVs fluctuation measured in an 
initial unbiased run may in fact turn to be suboptimal 
whenever the system is far from its initial configuration. 
In this respect the gradient-adapted choice is of great 
help. These methods simplify the choice of the input 
parameters for a metadynamics simulation, significantly 
reducing the time which is usually spent in trial and er- 
rors. Both approaches can be used with CVs which are 
bounded, cither by their intrinsic definition or by arti- 
ficially added potentials. In the latter case, the version 
based on time evolution of the CVs might be more ef- 
fective. Finally, they both can be straightforwardly com- 
bined with versions of MetaD which are based on multiple 
replicas 
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Appendix A: Bias divergence law 

The relation between the absolute free energy and the 
bias [Equation [2] and Equation 17 contains a shift C(t) 
which does not depend on s but grows with the time. As 
already noted, this constant is irrelevant in the calcula- 
tion of free-energy differences. However, it is instructive 
to compute explicitly its behavior in the long time limit. 
In the following we shall do it in the small-width ap- 
proximation, which is usually applied in the analysis of 
WTMetaD simulations. According to Equation [TJ] the 
time derivative of the bias is 



V(s,t) 



We -V{s,t)/AT yp^det^e- 



/ ds'e 



F(s>) + V(s> ,t) 



( A . 1} 

As in the long time limit the bias grows uniformly in s, its 
time derivative is equal to the one of C (t) . By combining 
Equation AT and Equation 15 one obtains 



C(t) 



tx e AT 



fds'{deta) s ; 



g T + AT 

By solving the differential equation Equation ~K2 



(A2) 



it can 



be seen that C(t) diverges logarithmically with time as 

C(t) 



lim , _ 

i^oo AT log t 



= 1 



(A3) 



We underline that this is just the limiting behavior of 
C(t). If its actual value is needed, for instance to align 
estimates of the free energy made at different times, the 
procedure illustrated in Ref. [3D] could be used to esti- 
mate it. Note however that in Ref. [3U] a different sign 
convention is used. 



Appendix B: Simulation details 

In all the simulations shown we employed alanine 
dipeptide (ACE-ALA-NME) in vacuum as a model 
system, with molecular interactions described by the 
CHARMM27jnij force field. This system (see Figure [TJ 
is a widely known benchmark for free-energy calculations 
[3"2Tl37| . Indeed, it displays two main basins, namely C7 eq 
and Cjn X separated by a sizable barrier of several ksT at 
300K. A timestep of 2 fs was employed and all the cova- 
lent bonds involving an hydrogen atom were constrained 
to the equilibrium distance by means of SHAKE [38 algo- 
rithm. A Langevin thermostat was used with a temper- 
ature of 300K and a damping factors of 5 ps _1 . NAMD 
2.8[55] molecular dynamics code was used and modified 
to add the adaptive shape WTMetaD capability. 



Appendix C: Error calculation 

Here we describe the procedure adopted to evaluate 
the error between two free-energy landscapes for the nu- 



9 



merical examples reported in the main text. 

We first note that WTMetaD for finite AT produces 
a more accurate histogram in low free-energy regions. 
Therefore we selected a reference region in CVs space 
defined by those points lying within a value of v free- 
energy units respect to the minimum of one of the free- 
energy surfaces, here termed reference free energy F r (s). 
In all the calculations performed the value of v was set 
to be 5 kcal mol -1 . 

In particular, being F(s) the free-energy surface whose 
error is required and s a point of the collective variables 
space, the error between the two surfaces is defined as: 



f s [F r (s) - F(s)]* 0(v-F r ( S ))ds 
V Is 0(v-F r (s))ds • 1 ' 

Here S is the multidimensional space in which the calcu- 
lation is performed, being $ and ^ for the Ramachandran 
plot, and 9 is a Heaviside step function. This equation 
amounts in calculating the average squared root differ- 
ence between the two free energies F r (s) and F(s) in the 
CVs space defined within v kcal mol -1 from the mini- 
mum in the reference free energy F r (s). The F r (s) is 
related to F r (s) by a rigid shift of the free-energy surface 
respect to the average value in the reference region: 



F r (s) = F r {s) 



f s F r (s)0(v - F r (s)) ds 
f s 6(v- F r (s))ds 



(C2) 



and a similar relation holds for the other free energy F(s) 
where, for consistency, the reference region is again de- 
fined on the reference free energy 



F(s) = F(s) 



fs F(s)6(v - F r (s)) ds 
J s 0(v- F r (s))ds ■ 



(C3) 



Appendix D: Supplementary Information: a simple 
illustrative toy model potential 

We want here to show how the need of using adaptive 
Gaussians arises from a simple one-dimensional Langevin 
model. Consider an energy landscape composed by 
the sum of three Gaussian potentials (a;o=-1.5, <7o=0-l, 
wo=-3, £1=1.5, o"i = 1.0, wi=-3, £2=0, (72=2.5, u>2=-3), 
namely: 



i=0 



exp 



(x Xi 

2af 



(Dl) 



and a particle of unitary mass that moves in this poten- 
tial with a temperature of 0.2 energy units. The timestep 



was set to 0.0025 time units and the friction was set to 
0.8 inverse time units. The particle is initially placed in 
the rightmost minimum (around 1.5) and the evolution is 
observed for 10 6 steps of Langevin dynamics. After this, 
the standard deviation of the position is calculated and 
used as width (<7tyT=0.26) for a Well- Tempered Meta- 
dynamics (WTMetaD) run during 2 x 10 6 steps of dy- 
namics. Additional parameters for the WTMetaD are a 
AT=2.8 and a energy deposition rate of 5 x 10 -5 energy 
units/step. In Figure [7] the result of this WTMetaD is 
shown along with the model potential (solid black line). 
A set of estimates of the underlying free-energy landscape 
(identical to the potential energy in this one-dimensional 
case) along the time of WTMetaD run are reported and 
aligned to the value of potential energy of the rightmost 
minimum. It is evident that the (Jwt which was suit- 
able for the larger minimum is not suitable for the min- 
imum on the left side thus seriously affecting both the 
free-energy difference between the minima as well as the 
barrier estimate. 
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FIG. 7: WTMetad run where awT=0.26 and stays constant 
along time. In solid black the energy profile of the model 
potential is reported. In color different estimates of the free- 
energy landscape along the WTMetaD run are reported. All 
the free-energy profiles are aligned to the free-energy value of 
the rightmost minimum in the model potential. 



On the contrary, in Figure[8]we report results obtained 
by setting the <jwt equal to the displacement of the 
system in the last 250 steps of the Langevin dynamics 
by means of our dynamical adapted algorithm (see main 
text). Here the narrow features of the free-energy are ad- 
equately resolved and the free-energy profiles and barrier 
result improved and present a faster convergence. 
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